% sampling from temporal spectrum
% creating the continuously random sequences of coefficients
function Gt = temporal_sampling(W_aj,params,C)
% W_aj: temporal spectrum of K-L function
% Output: Gt: time sequences of coefficients

coef_num = size(W_aj,3); % numbers of coefficents of Zernike poly

L = length(W_aj(:,:,1)); % 正频率长度 
N = (L-1)*2; % 数据长度 

fs = params.samplehz; % frequency of sampling
Gt = zeros(1,N,coef_num);
Xk = zeros(1,N);
for k = 1:coef_num
    % 《MATLAB数字信号处理的85个实用案例》---宋知用 P440
    Gx = W_aj(:,:,k);
    Gx(1,1) = 0; % 去除零频率分量
    Gxx = Gx;
    Gxx(2:L-1) = Gx(2:L-1)/2 * 2; % 由于本人编程中，在计算时间信号谱时并没有考虑负频率，所以计算出来的为完整的双边谱的正频率部分，再除以二会导致数值偏小了
    Ax = sqrt(Gxx*N*fs); % 1/t = fs/N; here: 1/t * N^2
    fik = pi*randn(1,L);
    Xk = Ax.*exp(1i*fik);
    Xk = [Xk, conj(Xk(L-1:-1:2))]; 
    Xm = ifft(Xk);
    Gt(:,:,k) = real(Xm)/sqrt(C(k,k));
end

% end of the function
end